1 Introduction 



Discrete kinetic theory, and notably the Lattice Boltzmann method has met 
with significant success in the recent years for the numerical simulation of a 
wide host of fluid phenomena, ranging from multiphase flows in grossly ir- 
regular geometries, up to highly turbulent homogeneous incompressible flows 
[|I|, [2|, [3]. Main drivers behind the method are ideal amenability to parallel 
computing, easy handling of grossly irregular geometries and physical sound- 
ness. 

The Lattice Boltzmann Equation (LBE) is more than a plain Navier- 
Stokes solver; actually, it is a minimal hyperbolic superset of the Navier-Stokes 
equations. In fact, in the process taking from the kinetic to the hydrodynamic 
level, along with standard primitive fluid fields such as density, speed and 
temperature, there is enough information left to describe the dynamics of 
a genuinely hydrodynamic stress tensor as well as that of a passive scalar. 
This property permits to extend the applicability of LBE techniques beyond 
'plain' hydrodynamics, with only minor modification to the basic discrete 
equation. To date, this feature has been exploited in various ways, to describe 
thermal convection, magnetohydrodynamics (0 and references therein) and 
also subgrid modeling ||. 

All of the above refers to the 24-speed FCHC (Face Centered HyperCube) 
lattice ||, the ancestor of modern Lattice Boltzmann schemes, introduced 
by D'Humieres, Frisch and Lallemand back in 1987. Recently, FCHC has 
been superseded by swifter Lattice BGK schemes 0, || which can boast a 
simpler (diagonal) collision operator as well as a significantly smaller number 
of discrete speeds (14 instead of 24). 

In spite of these undeniable advantages, we shall argue that the LBE- 
FCHC scheme still holds a certain interest for it provides the unique oppor- 
tunity to integrate passive and active scalars dynamics within the very same 
scheme tracking the flow variables, with no need to extend the corresponding 
set of discrete speeds. 

The scope of this paper is twodfold; first, to revisit the symmetry proper- 
ties of LBE-FCHC schemes, second to present a novel potential application 
of these symmetries to the field of reactive flow dynamics. 
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2 Reactive flow dynamics 



The set of hydro- chemical equations governing the dynamics of a S'-species 
reactive flow read as follows (k, I run over the spatial dimension): 

d t ps + di(p a ui) = diD s pdi(p s /p) + S s , s = l,S (1) 



d t pu k + di(pu k ui) = diupdiu k - d k p (2) 



d tP T + d k {pu k T) = d lX pd{T + Q (3) 

where p s is the density of the s — th species, p is the fluid density (p = 
J2 s Ps), T the fluid temperature, p the fluid pressure, u k the flow field, D, v 
and x t ne mass, momentum and energy kinetic diffusivity of the fluid. Here 
S s is the mass density input to species s from chemical reactions, and Q is the 
heat input rate as obtained by summing all over the reactions: Q = J2 r Qr^r 
, where ui r is the progress rate (moles/s) of reaction r. 

Combustion is governed by the competition between hydrodynamics (dif- 
fusion and convection) and chemistry; the relative strength of this competion 
is measured by a dimensionless parameter known as the Damkohler number, 
defined as Da = ^ where r h and r c are typical hydrodynamic and chemical 
timescales respectively, 

Low Damkohler 's characterize quasi-inert flows while high Damkholer's 
are associated with fast (r c << 77J chemical reactions. 

In this paper we shall be concerned with the 'fast' combustion limit, i.e. 
reactive flows in which chemistry proceeds much faster than hydrodynamics 
so that as soon as mixing is completed the fuel instantly burns ("mixed-is- 
burned"). 

In this limit, combustion tends to localize on thin regions ("reaction 
sheets" ) whose thickness becomes vanishingly small as the inverse Damkohler 
number goes to zero 0. This brings about a drastic simplification in the 
mathematical model because the focus can be shifted onto the topology of 
the reaction sheet. It can be shown that for the case of non-premixed flames, 
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i.e. flames where fuel and oxidizer are fed into the system in separate streams, 
and under the equidiffusional assumption (D = u = x), the topology of the 
reaction sheet is properly described by a passive scalar, say Z, evolving ac- 
cording to the usual advection-diffusion equation: 

pD t Z = d kP Dd k Z (4) 

where D t = d t + ud x + vd y + wd z is the fluid substantial derivative. 

The underlying idea is that while species transform one into another un- 
der reactive events, there are quantities (typically elemental mass fractions) 
which are unaffected by chemistry and serve therefore as slow modes en- 
slaving the system dynamics. As a result, in the limit where chemistry 
proceeds much faster than hydrodynamics, all species concentrations are 
'frozen-in' byproducts of the conserved scalars. With reference to a diffu- 
sive non-premixed flame, it is readily checked that a convenient choice for Z 
is provided by the fuel mixture fraction, defined as the mass ratio of fuel in 
the unburned mixture versus fuel mass ratio in the stream, 

Z = Y fu /Y fs (5) 

where Y stands for mass fraction and the subscripts u, s refer to the 
unburned and stream fuel respectively. 

For a single-step irreversible reaction of the form 
Fuel + Oxidizer — > Products 

combustion takes place around the surface where reactants and products 
come in stoichiometric proportion. The stoichiometric surface is defined by 
the relation: 

Z(x, y, z) = Z st = 1/(1 + AY f JY , s ) 

where A is the ratio of the fuel to oxidizer stoichiometric coefficients 



times the corresponding molecular weights [|K| . In the limit of infinitely fast 
chemistry, the stoichiometric surface represents a front of discontinuity for 
the spatial derivatives of the hydrodynamic variables. Within this limit, and 
with the additional assumption of steady state and adiabatic conditions, all 
quantities become a simple piecewise linear map of Z. With reference to lean 
(0 < Z < Z st , superscript "L") and rich (Z st < Z < 1), superscript "R") 
portions of the mixture, these mappings read as follows: 

Yf h {Z) = (6) 
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Y*{Z) = Y fs (Z - Z st )/(l - Z st ) (7) 

Y L b (Z)=Y 0S (l-Z/Z st ) (8) 

Y*(Z) = (9) 

T b L (Z) = T in + (Z/Z st )(T st - T in ) (10) 

T b R (Z) = T m + (±^-)(T st - T m ) (11) 

A — ^st 

where subscript " b" means " burned" . 



Here T in is the inlet stream temperature (we assume stream fuel and 
oxidizer at the same temperature) and the actual value of T st depends on the 
specific chemical reaction. Note that in the unburned portion of the mixture 
the dependence on Z disappears (T u = T in ). 

Within this approximation the heat source reduces to Q — Q r pAYf /Wf8(Z — 
Z st ) where Q r is the reaction heat (Joule/mole) and AYf represents the 
change of the fuel mass fraction between the unburned (Yf = ZYf s )) and 
the burned (Yf = Yf b ) state. The presence of the Dirac's delta reflects the 
instantaneous nature of the chemical reactions. 

3 Extended LBE for reactive flows 

The LBE is a minimal discrete Boltzmann equation reproducing Navier- 
Stokes hydrodynamics in the limit of small Knudsen numbers, i.e. particle 
mean free path much smaller than typical macroscopic variation scales [Q]. 
We shall refer to the 24-speed, single-energy, four- dimensional FCHC (Face 
Centered HyperCube) lattice defined by the condition 
E*= M 4 = 2 Ci k = (0,±1) 

The FCHC-LBE takes the form of a set of first-order finite difference 
equations: 

24 

fi(x k + c ik , t + 1) - fi(x k , t) = Y, Mfi - //) (12) 

i=i 
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where f,i = 1, 24 is a set of 24 populations moving along a corresponding 
set of discrete speeds and k = 1,4 runs over the spatial dimensions. 

The eq.(12) is naturally interpreted as a multi-relaxation scheme in which 
non-equilibrium gradients are brought back to the local equilibrium /? by 
scattering events mediated by the matrix A^. The local equilibrium popula- 
tions are chosen in the form of a discretized Maxwellian expanded to second 
order in the Mach number in order to retain convective effects. 

ft = + ZdkUk + 2(c ik cu - -8 k i)u k ui) (13) 

Under the constraint of point-wise mass and momentum conservations 
(J2iAij = J2iCikAij = 0) the equation (12) describes the motion of a four- 
dimensional fluid. The reason for working in four- dimensions is that no 
single-energy discrete lattice exists in three-dimensional space which ensures 
isotropy of the fourth order tensor T klmn = J2i CikCuC im c in . 

Three-dimensional motion is then recovered by quenching the propagation 
along the discrete speeds projecting upon the fourth dimension x 4 . Quench- 
ing the fourth dimension releases an excess of populations versus discrete 
speeds (24 against 18) introducing 'de facto' a degeneracy of the 'momentum 
eigenstates' (c^) propagating along the fourth dimension, since each of these 
links carries two populations with opposite speeds along £4. 

Generally, to the purpose of describing sheer hydrodynamics, there's no 
point of keeping two distinct populations ("doublet") along the same discrete 
speed. A 18 speed-18 populations (18sl8p for brevity) scheme is perfectly ad- 
equate for three-dimensional Navier-Stokes equations. If, on the other hand, 
doublets are retained, then the current J 4 is readily shown to be passively 
trasported by the three-dimensional fluid; this means that the 24sl8p model 
delivers the 3D Navier-Stokes equations plus a passive scalar. The same 
reasoning can be down- iterated to lower dimensions leading to 9s24p, 9sl8p 
(two-dimensional fluid plus one/two passive scalars respectively) and 3s24p, 
3sl8p 3s9p(one dimensional fluid and one/two/three passive scalars) in ID. 
A simple sketch illustrates the tree structure associated with dimensional 
quenching of the FCHC scheme: 
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24 4D: 4D fluid, passive scalar (24s24p) 



* 



18 6 3D: 3D fluid, passive scalar (18sl8p) 

* - 

9 9 6 2D: 2D fluid, passive scalar (9s9p) 

* - - - 2 (9s24p) 

3 6 9 6 ID: ID fluid, passive scalar (3s3p) 

1 (3s9p) 

2 (3sl8p) 

3 (3s24p) 

Since passive/active scalars are in heavy use in numerical combustion to 
describe fast chemistry regimes, we argue that the aforementioned properties 
point to the FCHC Lattice Boltzmann as to a potential candidate for the 
numerical description of combustion phenomena. 

In this paper, we shall focus on two-dimensional reactive flow dynamics 
with one passive scalar describing the fuel mixture fraction Z and a second, 
active, scalar accounting for temperature dynamics. This picks up the 9s24p 
model. It should be noted that within the "mixed-is-burned" assumption, 
the temperature field can be obtained directly from the knowledge of the 
mixture fraction Z via the piecewise linear mapping given by eq.(lO-ll). We 
choose however to track the temperature field independently of Z in order 
to gain a qualitative feeling for the validity of the adiabatic approximation. 

With these preparations, the reactive flow variables are defined as follows: 

P = Ei=i/i, pu = E?=iMi, pv = Ei=i/iQ2, pZ = Ei=i/iQ3, pT = 

Ei=l fi c U 

The second scalar T (temperature) is 'activated' by adding a source term, 
say q to all populations projecting upon the fourth dimension. A simple 
momentum budget fixes q as a function of the heat released Q. This is all 
we need to describe fast combustion. 



4 Numerical results 

As a test-case application, we shall consider a coflow two-dimensional methane- 
air laminar flame in a plane domain. The numerical set-up is as follows: 

Inlet: Coflowing Fuel (methane)-Oxidizer (air) at T = T in . The fuel 
enters the domain within an inner region centered about the midline and 
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width if/8. The mixture fraction at the inlet is Zi n . Outlet Zero cross- 
flow and zero axial derivatives of the streamwise velocity. Upper and lower 
boundary are kept at a constant wall temperature T w . Other simulation 
parameters are: inlet mixture fraction profile: Z in = 2.0, Z st G(y), G being a 
gaussian profile of width if/8, inlet temperature: T in = 0.0 (conventionally 
800 K), inlet flow: Poiseuille profile, maximum Ui n — 0.1, aspect ratio L/H = 
4. 

Our baseline computation is performed on a 256x64 grid and takes about 
1 [is I 'step/ 'site on a SUN Sparcl workstation, ninety percent of which due to 
hydrodynamics. Since no turbulence model has been implemented, the flow 
viscosity {y = 0.01) corresponds to a Reynolds number £ ~ ^ ne 
simulation is run over 3-5 axial residency times tr = 1.5L/Ui n = 3840. 

Chemistry is represented by the following one-step global reaction: 

CH A + 20 2 = 2H 2 + C0 2 ( 14) 

which is supposed to proceed instantly at the stoichiometric surface. 

For practical purposes, the reaction is assumed to occurr within a shell 
centered about Z = Z st = 0.057 with thickness Z t h = 0.2Z st . The heat 
source is assumed a (piecewise) constant within this shell. 

In Fig.l we show the equicontours of the mixture fraction at time t=4,000, 
corresponding to about one recirculation time. From this picture the typical 



shape of a laminar diffusion flame is apparent JTT], [12[ . Note that the sto- 
ichiometric surface is confined to the initial region of the domain (x ~ 50 
lattice units). For a given set of physical-geometrical parameters, this loca- 
tion, namely the flame length Lf, depends solely on Z in , i.e the degree of 
mixture richness. 

In Fig 2 we show the mixture fraction at the midline y = H/2 as a function 
of the streamwise coordinate x at two different time snapshots t = 2000, 4000. 
The advancement of the mixture front is well visible along with its diffusive 
spreading. As expected, after an initial diffusive transient lasting about 
Tr/ Re ~ 100 time steps, the tip of the front proceeds on a quasi-convective 
regime at speed Uf = 2C/j n /3. The flame however stabilizes around the 
stoichiometric surface long before the tip comes to the outlet of the domain 
because L f is much shorter than the longitudinal span L of the computational 
box. 

Within the equilibrium chemistry assumption, mid-plane reactions (y = 
H/2) should peak around the stoichiometric line x = x st , where Z(x st ,y = 
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H/2) = Z si . In our simulation, x st — 52 as highlighted by the horizontal 
line in Fig. 2. This is the region where a substantial temperature buildup is 
expected to occurr. Such a prediction is indeed reflected by Fig 3, where the 
median temperature profile is reported as a function of the axial coordinate x. 
In line with the behaviour of the mixture fraction, the temperature field shows 
a peak centered about the stoichiometric region, with a top temperature of 
about 2100 K, slightly below the stoichiometric temperature T st = 22Q3K. 
For the sake of comparison, the adiabatic profile deduced by the theoretical 
"frozen-in" mapping T(Z) computed with Z at t = 4000 is reported right 
above (dotted line on top). From this figure, we see that a relatively good 
match is obtained in in the post-flame region, while in the pre-flame region 
the adiabatic profile provides a gross overestimate of the temperature field. 
This is not surprising, since the near-inlet region is characterized by strong 
gradients which make the adiabatic assumption rather questionable |T3[ . 

Finite-rate chemistry effects are typically highlighted by means of temperature- 
mixture fraction scattergrams, such as the one presented in Fig. 4. This 
figure, referring to the transversal section x = 40 after 4, 000 time steps, 
conveys a clear idea of the cooling effect produced by diffusive processes. 

Part of these diffusive processes is likely to result from lattice discreteness 
as well, and in particular from the finite-thickness Z t h of the reaction region. 
These numerical effects can be controlled by a careful tuning of the discretized 
version of the Dirac's delta function 8(Z — Z st ), but no systematic studies 
have been conducted in this respect. 

In Fig 5 we present the transverse temperature profile at three different 
axial locations x = 1 (inlet), x = Lf/2 (mid-flame) and x = 128 (mid- 
domain). As expected, a doubly-humped profile is observed within the flame 
extension, with the two peaks centered about the stoichiometric surface (com- 
pare with Figure 1). Beyond the flame front a typical gaussian profile is re- 
covered. The overall structure of the temperature field is finally represented 
in Fig 6, where the equi contours of the temperature field are displayed. From 
this figure, the transition between double to single humped transverse profiles 
is even more apparent. 

In conclusion we have shown that minor modifications to the basic FCHC 
lattice Boltzmann scheme permit to deal with reactive flows in the limit of 
infinitely fast chemistry. Given the fact that chemistry is intrinsically local in 
space and Lattice Boltzmann is a proved parallel performer on its own even 
for inert flows, excellent performances on parallel machines may be easily 
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anticipated. This said, it is worth emphasizing that this work is still very 
preliminary in character and several extensions are required before quan- 
titative comparisons with state-of-the art numerical combustion techniques 
can be attempted. Among others, important future developments are the 
following ones: 

1. Two-dimensional laminar to three-dimensional turbulent flames 

2. Infinitely fast to finite-rate chemistry 

3. Nearly incompressible to compressible flows 

The first item does not require any modification to the scheme presented 
here, but just higher numerical resolution. To convey a quantitative idea 
of how far we can go with direct simulations, let us mention that a purely 
hydrodynamic model (18sl8p) for isothermal turbulent channel flows is cur- 
rently achieving Reynolds numbers of about 3, 000 on the massively parallel 
computer Quadrics, with a sustained speed over 10 Gflops ||I4[. The nu- 
merical exploration of this regime should provide valuable informations on 
the morphology of turbulent flames. Higher Reynolds numbers will call for 
suitable turbulence models. These can easily be incorporated within the 
LBE method so long as simple algebraic turbomodels such as Smagorinski, 
Baldwin-Lomax and Renormalization Group models are used. All what is 
needed is to let the scattering matrix Aij a local function of the stress tensor, 
i.e. a pointwise function of the populations ||. 

Moving beyond the fast chemistry regime requires the development of 
a suitable multispecies transport scheme able to handle several scalars in 
three dimensions. Since in the conventional Lattice Boltzmann formalism 
each extra-scalar requires six moving populations, a naive extension of the 
present scheme does not look very appealing. However it seems possible to 
devise modified LBE schemes which are able to track passive scalars using 



only population per species ||18j| . An alternative option is to couple the 
present LBE scheme to 'conventional' grid-based methods. 

At variance with the previous two, the third item does involve substantial 
changes to the scheme proposed in this work. 

The scheme presented in this paper belongs to the simplest class of com- 
bustion models; those in which density is not able to keep up dynamically 
with significant temperature excursions. So, in principle, the method should 
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only be applied to 'cold' flames with little heat release in which the main 
focus is on flame morphology (i.e. the scenario pertaining to item 1. in the 
aforementioned list of developments). 

The next class of problems to be adressed is low-Mach reacting flows in 
which density is allowed to respond to temperature changes over a signifi- 
cant dynamical range of values. The possibility to model this class of flows 
by means of appropriate phenomenological extensions of the present LBE 
scheme is currently under investigation. 

Finally, one should consider the class of fully-compressible, high-Mach 
flows in which density is supposed to respond to both temperature and pres- 
sure changes. We believe such class is definitely out of reach for the present 
LBE model. 

In fact, in order to capture full thermo-hydrodynamic behaviour the dis- 
crete speeds must necessarily accomodate multienergy levels. A few models 
in this class have already been proposed |TJ|, |I6 |, but they still need to be 



perfected in many respects. First, they involve many discrete speeds (40 
in 3D), let alone the additional degrees of freedom explicitly required to 
incorporate passive/active scalar dynamics; second, and more importantly, 



it appears like their stability is still matter of concern [I7|. Further re- 
search work is definitely required to clarify whether discrete speed models 
(of reasonable complexity) will ever be able to perform competitively for full 
thermo-hydrodynamic combustion applications. 

Meanwhile, we believe that the model presented here, although clearly 
limited in scope, can provide a useful starting point to explore phenomeno- 
logical alternatives to fully-fledged discrete speed models for the numerical 
simulation of (some class of) combustion phenomena. 
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Figure captions 



Fig.l: Isocontours of the mixture fraction field Z(x, y) after 4, 000 time 
steps. The stoichiometric line corresponds to Z — 0.057. 

Fig. 2: The axial profile of the mixture fraction at the midline y = 128, 
after 2, 000 and 4, 000 time steps. The solid line is associated with the 
stoichiometric value Z st = 0.057. 

Fig. 3: The axial profile of the temperature field at the midline y = 128, 
after 2, 000 and 4, 000 time steps. The solid line is the initial profile and 
the curve labeled 'Tadiabatic' represents the adiabatic profile computed 
with the mixture field Z at t — 4, 000. 

Fig. 4: Temperature versus mixture fraction scatterplot at x = 40 at 
time t = 4,000.. 

Fig. 5: The transverse temperature profile at three axial locations, x — 1 
(inlet), x = 25 (mid-flame) and x = 128 (mid-domain) at t — 4, 000.. 

Fig. 6: Isocontours of the temperature field after 4, 000 time steps. 
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